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Abstract 

The definition of artificial immunity, realized through vaccinations, is nowadays a practice widely developed in 
order to eliminate cancer disease. The present paper deals with an improved version of a mathematical model 
recently analyzed and related to the competition between immune system cells and mammary carcinoma cells 
under the action of a vaccine (Triplex). The model describes in detail both the humoral and cellular response of 
the immune system to the tumor associate antigen and the recognition process between B cells, T cells and 
antigen presenting cells. The control of the tumor cells growth occurs through the definition of different vaccine 
protocols. The performed numerical simulations of the model are in agreement with in vivo experiments on 
transgenic mice. 



Background 

The immune system (IS) is a complex system of organs, 
cells and molecules whose main task is to protect living 
organisms from external pathogens such as viruses and 
bacteria. Nevertheless the effectiveness of the IS in tumors 
disease is nowadays under discussion among biologists 
and physicians. As stated by the immunosurveillance the- 
ory [1,2], biotechnology engineered naked mice (mice 
without immune system) show the developing of multiple 
variants of malignant tumors that are not usually visible in 
wild mice, thus suggesting that the immune system plays 
an important role also against tumors. Indeed most 
mutated malignant cells are recognized and eliminated by 
immune system mechanisms right after their birth, and 
tumors that usually arise are indeed poorly immunogenic 
tumors, originating from malignant cells which escape 
from immune surveillance. Some tumors are caused by 
exogenous factors (e.g., smoke for lung cancer), and the 
elimination of the exogenous cause would in theory pre- 
vent the risk of developing the tumor. However many 
other tumors are caused by endogenous factors and their 
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developing cannot be easily predicted and controlled. 
Among human cancers, the mammary carcinoma repre- 
sents a major cause of concerns in women, since it belongs 
to the class of endogenous cancers which escape immuno- 
surveillance of the IS. 

The risk of appearance of mammary carcinoma is 
usually estimated by analyzing the family history of can- 
cer, and breast cancer screening in young women is 
highly recommended since the achievement of earlier 
diagnosis could greatly improve the outcomes of the 
treatment. Strong family history of cancer usually entitles 
higher risks of developing the tumor, thus suggesting 
that tumor hereditary is encoded into the DNA. Some 
gene tests such as the genetic screening for the BRCA 
genes [3] are nowadays possible and may determine the 
risk of cancer. Indeed the analysis of the genome of indi- 
viduals will be useful to better estimate the risk of cancer. 

Biologists and physicians are exploring novel immuno- 
preventive treatments that can avoid the development of 
breast cancer in patients with high risks of malignant cell 
mutations. Among others, Lollini et al. [4] have devel- 
oped a cellular vaccine, called Triplex, which is able to 
elicit complete prevention of mammary carcinogenesis in 
HER-2/neu transgenic mice. Triplex combines three dif- 
ferent elements (the tumor antigen with two adjuvants) 
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that stimulate the immune system response with different 
actions [4] . Vaccine cells have been engineered to present 
and release the tumor associated antigen p-185 (that is 
also the oncogene of the tumor) with the addition of 
Allogeneic MHC-class I molecules to easier recognizing 
by cytotoxic T cells. Moreover, thanks to transduction of 
interleukin-12 genes, they release interleukin-12 mole- 
cules that have a broad range of costimulatory functions 
in boosting the immune response against tumors. 

It is worth stressing that differently from the vaccine for 
virus or bacteria, cancer vaccines require repeated admin- 
istration for the the entire life of the host. This is due to 
the low antigenicity of the cancer cells, the capability to 
escape the immune system surveillance. Moreover present 
cancer immunoprevention research is unable to find better 
vaccines able to assure complete, long-term protection. 

The repeated administration of the vaccine, realized 
with the aim to increase the antigenicity of the tumor 
associated antigens, maintains the immune system 
response ready against newborn cancer cells. However, 
even if vaccines are usually less toxic than standard 
drugs, uncontrolled administration of the vaccine can 
induce undesirable effects such as autoimmune diseases. 
Therefore the optimization of the vaccination protocol 
constitutes a fundamental and open problem. 

In the in vivo experiments it is not usually possible to 
reach an optimum vaccination protocol that maximizes 
the efficacy of the tumor depletion on the one hand and 
minimizes the risk of side effects on the other hand, 
because of the large variability cases. Indeed vaccination 
protocols are usually determined heuristically basing on 
best practice and previous experience. Moreover the cost 
of in vivo experiments can be prohibitive. 

In order to understand whether it was possible to gain 
complete prevention of mammary carcinogenesis with 
fewer injections, a (multi) agent-based model named Sim- 
Triplex [5] has been developed. It is worth noticing that 
SimTriplex has been also employed for other pathologies 
[6-10]. However agent-based models do not allow the 
development of asymptotic analysis of the competition 
and an easy investigation of parameters' space. 

Different mathematical tools have been developed in 
order to model complex biological systems and among 
others, immune system-cancer competition. The most 
famous approach is the ODE-based model where the over- 
all system is decomposed in different cell populations 
whose time evolution is depicted by solutions of a non- 
linear ODE system (nonlinear terms take care of the inter- 
actions among two or more cell populations), see paper 
[11] for a review of ODE models available in the literature 
and [12] for a comparison between ODE models with and 
without delay. 

Kinetic theory models have been also proposed for the 
immune system-cancer competition. These models 



consist in partial integro-differential equations and allow 
both the modeling of proliferative/destructive events 
and the modeling of mutations occurring in the onset of 
tumor, [13]. Further modeling approaches for the 
immune system-cancer competition include cellular 
automata, agent-based models, see the recent expository 
paper [14]. 

Most of the mathematical models of the IS summarize 
the response of the immune system in a single popula- 
tion of cells, named effector cells, which perform the 
task of destroying cancer cells. This simplifying assump- 
tion allows to reduce the complexity of the dynamics of 
immune system but it neglects the recognition process 
that occurs among the different cells that constitutes the 
response of the IS to the tumor antigen. 

The ODE-based model proposed in this paper has been 
derived from a biological conceptual model that is a good 
representation of the biological scenario (see Figure 1). 
The model takes into account both the humoral and cel- 
lular response of the immune system and the recognition 
process that involves the following entities: vaccine cells 
(VC), cancer cells (CC), tumor associated antigens 
(TAA), Plasma B cells (B), thymus cytotoxic lymphocytes 
(TC), thymus helper (TH) lymphocytes, antibodies (AB), 
interleukins 2 and 12 (IL2 and IL12), and antigen pre- 
senting cells (APC). A simplified version of the mathema- 
tical model proposed in the present paper has been 
analytically investigated in [15,16]. The simplified model 
does not include the role of the associated antigens, 
plasma B cells, interleukins 2 and 12 and the antigen pre- 
senting cells. Therefore the mathematical model of the 
present paper is a robust extension, from the biological 
viewpoint, of the model analyzed in [16]. In the present 
paper we restrict our attention to the comparison of 
numerical solutions of the model with in vivo experi- 
ments and sensitivity analysis of the model parameters. 
The model described in this paper can be also developed 
in order to take into account many biological phenom- 
ena, like chemotaxis, spatial cell dynamics or cluster for- 
mation. If spatial cells dynamics needs to be included, 
one can use the mathematical framework of the kinetic 
theory for active particles, see [17] and the reference 
therein. According to the latter framework, cells are 
grouped in functional subsystems which express a speci- 
fic strategy (called activity) and the time evolution of the 
subsystem is represented by a distribution function over 
the cells microscopic state (position, velocity and activ- 
ity). In this framework the mathematical modeling of the 
chemotaxis phenomenon and the formation of tumor at 
tissue scale can be included as shown in [18] and [19,20]. 

The present paper is organized as follows: Section "The 
Triplex vaccine in vivo experiments" briefly deals with the 
phenomenological analysis of the biological system. Sec- 
tion "The ODE-based model" is devoted to the description 
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VC are killed by 



AB 




Figure 1 Conceptual model of the in vivo experiment. On top vaccine cells (VC) are administered through intravenous injections, and then 
recognized by Cytotoxic T cells (TC) and Antibodies (AB) that kill them. Killed VC release both lnterleukin-1 2 (IL12) and Tumor associated 
antigens (TAA). TAA are captured by antigen presenting cells (APC) and then presented to T helper cells (TH). IL12 stimulates both TH and TC 
actions. TH release interleukin-2 (IL2) which boosts TH, TC, and B actions and stimulates B cells to differentiate into plasma B cells (B). B release 
AB, and both AB and stimulated TC kill cancer cells (CC), which further release TAA. 



of the ODE-based model. Section "Sensitivity analysis" 
introduces the sensitivity analysis technique. Section 
"Results and discussion" compares, by means of numerical 
simulations, the mathematical model with the computa- 
tional model SimTriplex. Finally Section "Conclusions" 
concludes the paper with a critical analysis and research 
perspective of the model. For interested readers Additional 
File 1 presents a simplified version of the model by cou- 
pling the differential model with an algebraic model. 

Materials and methods 

The Triplex vaccine in vivo experiments 

This section briefly deals with the in vivo experiment car- 
ried on BALB-neuT neu virgin female mice groups which 



over-express the activated rat HER-2/neu oncogene. The 
description does not pretend to be exhaustive from the 
biological point of view but highlights the essentials of 
the experiments in order to motivate our study. 

The Triplex vaccine has been obtained from a mam- 
mary carcinoma of a FVBneuN #202 (H-2 9 ) mouse, 
transgenic for the rat protooncogene c-neu, and com- 
bines different stimuli: 

♦ The pl85neu oncoantigen; 

. The H-2 q MHC molecules (allogeneic for H-2 d 
BALBneuT mice); 

• The interleukin-12 (vaccine cells are engineered 
with the genes coding for murine IL-12). 
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The experiment starts at the sixth week of age, where 
BALB-neuT mice start the vaccination protocol. Mice 
are divided in different groups, one for control 
untreated group, and one for each protocol tested. All 
vaccine protocols that have been tested are built upon 
the same 4-week cycle which consists in twice-weekly 
intra peritoneum vaccinations (Tuesday and Friday) for 
the first 2 weeks followed by 2 weeks of rest. 

The Prophylactic, lifelong Chronic vaccination proto- 
col of cancer-prone HER-2/neu transgenic mice with 
cells expressing HER-2, allogeneic MHC antigens and 
IL-12 demonstrated able to completely prevent the 
onset of mammary carcinoma. The Early vaccination 
protocol (which counts only three 4-week cycles at the 
beginning of the experiment) produces a significant 
delay in the onset of tumors, but all mice eventually 
succumbe to mammary carcinoma. Other tested proto- 
cols demonstrated much less effective, with little or no 
gain in efficacy when compared to untreated control 
mice. 

It is worth stressing that maximal prevention against 
mammary carcinoma required all the three vaccine com- 
ponents (HER-2/neu, allogeneic MHC antigens, and 
IL-12) and was due to the induction of both cellular and 
humoral immune responses. Although cellular and 
humoral immune responses are taken into account in the 
vaccine administration, the relative importance of antibody 



subclasses for successful cancer prevention indicates that 
humoral immune responses is more important than cellu- 
lar responses driven by cytotoxic T cells [4] . 

Recent investigations [21] show that the Triplex vaccine 
progressively looses its efficacy with the advancement of 
tumor progression, both in terms of tumor incidence and 
multiplicity (see Figure 2). In particular, tumor develop- 
ment is remarkably delayed in mice receiving the early 
protocol with respect to untreated mice, whereas protocols 
started later have produced only a negligible delay. 
Furthermore in vivo tests show that the Triplex vaccine is 
ineffective against larger tumor targets. Thus, the triplex 
vaccine demonstrates very effective at preventing mam- 
mary carcinoma onset in tumor-free mice but is ineffective 
against established local tumors. 

It should be therefore clear that any vaccination proto- 
col should be started early enough to avoid carcinoma in 
situ formation. On the other hand it should be advisable 
to minimize the number of administrations in order to 
both maintain complete efficacy and reduce the risk of 
any undesirable effect. In order to help biologists in find- 
ing better vaccination protocols, a (multi) agent model 
named SimTriplex has been developed in [5]. The model 
has been inspired by the work of Celada and Seiden [22] 
and uses an approach that models ab initio the interac- 
tion and diffusion kinetics of each relevant biological 
entity. SimTriplex has been tuned with the in vivo 




Tumor Progression (time) 

Figure 2 Triplex vaccine efficacy measured in in vivo experiments with respect to the advancement of the tumor. The abscissa 
represents the main temporal stages of tumor progression: from atypical hyperplasia up to mammary carcinoma. The ordinate shows the rate of 
inhibition of tumor burden entitled with the use of the vaccine. The red line represents the achievable efficacy of the Triplex vaccine in 
preventing the tumor burden whether the first protocol administration is delayed at successive stages of tumor progression. 
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experiments and demonstrated able to coherently repro- 
duce the behaviors of the entities involved in the in vivo 
immunoprevention experiment. In addition the use of 
SimTriplex as a predictive tool yielded to encouraging 
results [23]. 

The ODE-based model 

The competition between immune system cells and can- 
cer cells reminds the well known Predator-Prey (PP) 
model described by Lotka-Volterra equations. There is a 
population of prey, represented by the cancer cells, with 
an infinite set of food resources (nutrients coming from 
the host blood) and, differently from the classical PP, 
multiple populations of predators (the effectors cells) 
cooperate through cell-cell and cell-molecule interactions 
to neutralize the prey. Differently from the classical pre- 
dator-prey models, predator survival does not depends 
on the number of prey, since predator populations exist 
normally and, in absence of the prey, their number oscil- 
late around given equilibrium levels. 

If cancer cells are able to escape immunosurveillance, 
the cancer takes over, tends to compete with the healthy 
cells for nutrients, and could be able to kill the host. How- 
ever the immune system response can be helped in recog- 
nizing harmful cells by an external agent represented in 
our case by a vaccine. The induced immune response is 
the result of a complex network of interactions between IS 
cells which mainly depends from cells receptors. IS cells 
which present, through their receptors, specific tumor 
antigens can trigger a complex process whose final result 
is the eradication of the tumor. Specific interactions, 
which involve cell receptors, cannot be described with an 
ODE population model, see [5], however if we assume that 
the vaccine cells activate IS cells at given ratios we can 
model the subsequent immune response of activated cells. 

The network of organs, cells and molecules involved 
in the immune system is very large. In the model of the 
present paper we include only the entities recognized as 
fundamental for the biological process. We also assume 
that the IS-cancer competition occurs only in one 
hybrid organ which includes all the physical involved 
compartments (peritoneum, mammary gland, lymph 
nodes and so on). 

Both the humoral and cellular responses of the immune 
system are taken into account with their main entities, 
plasma B cells (B), cytotoxic and helper T lymphocytes 
(TC and TH), antibodies (AB), antigen presenting cells 
(APC) and interleukins 2 and 12 (IL2 and IL12), while the 
pathogens are represented by the cancer cells (CC), the 
vaccine cells (VC) and the P-185 tumor associated anti- 
gens (TAA). The interactions among the various entities 
of the immune system network, the external stimulus of 
the vaccine and the cancer cells are modelled by a system 
of ten nonlinear ordinary differential equations whose 



variables are summarized in table 1; In Figure 1 we show 
the conceptual model for the biological problem. 
Model parameters 

The model contains 44 parameters which have a specific 
biological meaning. These parameters, assumed as con- 
stants, modify the rate of variations of the populations 
due to natural death, interactions with other population 
and release of new quantities. Accordingly: parameters 
referring to natural death of entities are denoted by 
where i identifies the population under consideration; 
parameters referring to the interaction between popula- 
tions i and / are identified by a, finally parameters 
referring to releasing processes are identified by y ; ; , 
where i refers to the released entity and /' to the releas- 
ing entity. 
Vaccine cells 

The vaccine cells dynamics is described by equation (1). 
Vaccine cells {VC) are injected into the host through 
intraperitoneal vaccination with a predefined dosage. 
The inoculation of the vaccine cells is modeled by a 
function k in (t, q) which adds q vaccine cells to the cells 
in the host at time t whether a that time an injection 
was scheduled. Since vaccine cells come from the exter- 
nal, this term represents the only source element in the 
equation. Vaccine cells die for multiple causes, such as 
natural death (term -^VC), are inhibited by cytotoxic T 
cells that recognize vaccine cells thanks to their allo- 
genic-MHC class II molecules (term -a 19 TCVC), or by 
specific antibodies that are able to directly kill vaccine 
cells by complement mechanism (term -a 17 ABVC). 

— — = fej„(t, cj) - fiiVC - {a 19 TC + a 17 AB)VC (1) 
at 

P-185 tumor associated antigens 

Equation(2) models tumor associated antigens dynamics. 
Tumor associated antigens (TAA) can be released by 
dead or killed vaccine or cancer cells. Accordingly we 
suppose that the number of released antigens is propor- 
tional to both the number of vaccine cells and the 



Table 1 Model variables 


Variable 


Description 


Short name 


*i 


Number of injected vaccine cells 


VC 


Xl 


Number of P-185 tumor associated antigens 


TAA 


*3 


Number of activated B cells 


B 


X 4 


Number of activated T helper cells 


TH 


*5 


Number of interleukin 12 molecules 


I LI 2 


X6 


Number of interleukin 2 molecules 


IL2 


x? 


Number of released antibodies 


AB 




Number of cancer cells 


CC 


X 9 


Number of activated cytotoxic cells 


TC 


x 0 


Number of activated antigen presenting cells 


APC 
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number of cancer cells (VC and TAA respectively) that 
are inhibited (terms y 2 i(a 19 TC+ a^AB+^VC and y 28 
(a 8S +a &9 TC)VC). These two terms represent the source 
elements of the equation. Antigens are subjected to nat- 
ural degradation (term (i 2 TAA) and phagocytosis by 
antigen presenting cells (-a 2 oAPCTAA), such as dendri- 
tic cells, macrophages and B cells. Moreover antibodies 
can bind to free antigens producing immune complexes 
(term -a 27 ABTAA). 

dTAA 

= Y2\{a a TC + anAB + m)VC + Y2&( a as + a m TC + a sl AB)CC+ q\ 
- (lJ.2 + a 20 APC + tx 17 AB)TAA 

T helper cells 

The key role of T helper cells is to stimulate both the 
humoral and cellular branches of the immune response 
by direct receptor binding or the releasing of specific 
cytokines that boost the immune response, such as 
interleukin 2. T helper cells are activated by specialized 
APC such as dendritic cells, macrophages or presenting 
B cells which present major histocompatibility class II/ 
peptide complexes. 

Since we are dealing with a monoclonal model, pre- 
sentation is not directly modeled, so the percentage of 
activated T helper cells is estimated from the number of 
antigen presenting cells (APC) existing in the system 
(term y 40 APQ. Interleukins 12 and 2 (/LI 2 and ZL2) 
contribute to stimulate T helper cells priming and dupli- 
cation (termsa 46 (nH7)THanda45(7ni^7)TH). It is 

worth noting that, since interleukin 2 is released by T 
helper cells, these cells are able to self-stimulate their 
activities. The death factor is modeled by -ft 4 TH term. 

dTH / IL2 \ I 1L\2 \ ,„. 

Plasma B cells 

Plasma B lymphocytes may absolve to multiple functions 
in building the immune response chain against patho- 
gens. In a first stage they can act as specialized antigen 
presenting cells, by recognizing pathogens through their 
specialized "Y-shaped" receptors, and can then present 
peptidic sequences to T helper cells. As a consequence 
of a successful interaction with T helper cells, they can 
be stimulated to differentiate into plasma B cells, which 
release antibodies with the same receptors shape, or B 
memory cells, which readily act against new appearance 
of previously encountered pathogens. Since there is no 
in vivo experimental evidence of B memory cells, as also 
suggested by the need of a chronic vaccination to 
achieve complete protection against tumor onset, we 
decided to do not include for now memory B cells into 
the model [4]. 

In equation(4) we only consider the behavior of the B 
cells population (B) that has been activated by T helper 



cells (TH) positive feedback (term y 34 TH) and is there- 
fore able to release specific antibodies against cancer 
cells. We include the B as APC function in equation 
(10). Interleukin 2 (IL2) released by T helper cells plays 
an adjuvant role in stimulating B cells duplication 

(term «36( j/"^ )B) . Death is modeled by -p 3 B term. 
dB ( IL2 \ 

^-^ TH + ^{iL2-77j B -^ B (4) 

Interleukin 12 

Interleukin 12 (/L12) is mainly introduced through vaccine 
administrations, so it depends on the vaccine dosage. In 
previous in vivo experiments [24] interleukin 12 was intro- 
duced separately, but after transduction of IL2 genes 
inside vaccine cells [4], it is released by killed vaccine cells, 
so it is proportional to the number of killed vaccine cells 
(term y 51 (a 19 TC + a 17 AB + fi^VC). IL2 is subjected to 
normal degradation (-n^LYl) and it is partially absorbed 
for mitotic and stimulation signals by cytotoxic and helper 
T cells priming (terms -a 59 TCIL12 and -a 54 THIL!2). 

—Jf- = X5i(ai9TC + a„AB + /Xi)VC - (a 54 TH + a 5 ,TC + /t 5 )JL12 (5) 

Interleukin 2 

Interleukin 2 is mainly released by T helper cells (term 
y 64 TH). As previously stated, interleukin 2 stimulates T 
helper priming, and primed T helper cells produce 
further interleukin 2. It is subjected to normal degrada- 
tion (-fi 6 IL2) and it is partially absorbed for mitotic and 
stimulation signals in cytotoxic T cells priming (term 
-a 69 TCIL2) and B cells duplication (term -a 63 BIL2). 

dIL2 

— — = y 64 TH - {a 63 B + a 69 TC)IL2 - \± 6 \L2 (6) 
dt 

Antibodies 

Antibodies represent the main result of the humoral 
immune response. Antibodies (AB) are released by plasma 
B (B) cells (term y 13 B) and are subjected to normal degra- 
dation (modeled by -ft 7 AB term). More-over they disap- 
pear in absolving their functions: binding to specific 
targets, i.e. antigens (term a 72 TAAAB), cancer and vaccine 
cells (terms a 7S CCAB and a 7l VCAB, respectively). 

— — = y 73 B - [a 7& CC + a n VC + a 72 TAA]AB - /i 7 AB (7) 
at 

Cancer cells 

Cancer cells growth (CC) is modeled through the term 

[(i — k — assj CC . The term -a 88 CC is used to 

take into account CC killing by other immune system cells 
that are considered of minor importance for the process 
and are consequently not explicitly modeled, such as 
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Natural Killer cells which can kill cancer cells that under- 
express the major histocompatibility class I complex. The 
term p models the continuous production of newborn 
cancer cells. Due to transgenic nature of HER-2/neu mice 
new cancer cells are in fact continuously introduced into 
the host. The other terms describe cancer cells death 
mainly due to antibodies (term -a &7 ABCC) and cytotoxic 
T cells (term -a &9 TCCQ actions. As matter of a fact acti- 
vated cytotoxic T cells can kill cancer cells by direct cyto- 
toxicity and specific immunoglobulins can kill cancer cells 
by complement and other mechanisms. 



dec iy cc\, 

— — = 1 \k - G-88 

at LV CmaxJ 



CC - (a 89 TC + a s7 AB)CC + p (8) 



Cytotoxic T cells 

Cytotoxic T cells priming {TC) depends mainly on vac- 
cine cells {VC). Vaccine cells are engineered with allo- 
geneic major histocompatibility class I complex in order 
to easier presentation (term y 91 VC). Duplication is 
instead indirectly stimulated by T helper cells through 

the release of interleukin 2 (term a 96 ( IL 2+ S96 )TC)). Nat- 
ural death is modeled with the term -fi 9 TC. 



dTC 
dt 



y 91 VC + Gf 9 6 



(—) 

\IL2 + s 96 J 



TC — flgTC 



(9) 



Antigen presenting cells 

With the term antigen presenting cells we indicate a class 
of different types of cells, such as dendritic cells, macro- 
phages, but also B cells, whose focal mission is to recognize, 
capture, and process antigens in order to present small 
antigenic sequences named peptides in conjunction with 
MHC class molecules to both cytotoxic and helper T cells. 

Antigen Presenting cells (APC) are then depending on 
the quantity of the antigens that have been released 
(term y 02 TAA), and can die (term -ftoAPQ. 



dAPC 
dt 



y 02 TAA - ixqAPC 



(10) 



We thus designed the following set of ten non linear 
ODEs that is able to model the considered system of 
cell populations and interactions: 



iiVC 
d'lAA 



da _ 
IT ~ 
dm 
dt 

dILll 



dAii 

dt ' 
dCC 

dt 
dTC 

"3T : 

d\l\ 
dt 



-- hn(t, q) - IM VC - (ai 9 TC + a 17 AB) VC 

= Y2i(aiiTC + a 17 AB + /j.i)VC + x 28 (£«88 + a&gTC + a s7 AB)CC+ 

-(^ 2 + a 20 APC + a 27 AB)TAA 
y 34 TH + a 36 (jJ^) B - n 3 B 

- y i0 APC + a 46 (n§^) TH + a 45 (^) TH - , M TH 
= Y5i[oci9TC + a\ 7 AB + n\)VC — (a^TH + as^TC + fi5)lL12 

= y^TH — (a^B + a^TC)IL2 — ^i^lLl 

■- y 73 B - [a 7s CC + a vl VC + a 72 TAA]AB - fi 7 AB 

-- [(l - ^) k - a 88 ] CC - (a 89 TC + a 87 AB) CC + p 
Ys\VC + a 96 (jnnlr) TC - M9 TC 

= Y02TAA - n 0 APC 



Since we consider populations that are activated by 
vaccine administrations, we set the following initial con- 
ditions: 

VC(0) .IM(0) -B(0) .'J'll(O) -1112(0) -;t2(0) -AB(0) -CC(0) -TC(0) -A1"C(0) - 0 

The parameters in the model have been derived from 
literature, from measurements made during the in vivo 
experiment and from the SimTriplex model. Some para- 
meters which belong to the class of free parameters that 
any model has, were chosen into reasonable rages in such 
a way that the model was able to reproduce in vivo mean 
survivals for the untreated, early, and chronic vaccina- 
tions using a trial and error technique with mean-square 
evaluation, see table 2. 

During the in vivo experiment, biological dynamics is 
observed in time slices that are not smaller than eight 
hours. For this reason, we set the simulation time step 
equal to (A(t) = 8 hrs). This biological motivation also 
determined the SimTriplex time-step. The choice of the 
physical time-step allows to compare the results of the 
two models. Both models are supposed to simulate the 
dynamics of entities inside a volume of lfil, which corre- 
sponds to a small portion of mammary gland of mice. 

Sensitivity analysis 

In order to understand which parameter may be consid- 
ered fundamental in this process, it is significant to inves- 
tigate the sensitivity of the model to the alteration of the 
parameters. Choosing a parameter in a suitable range 
while retaining fixed the others, represents the classical 
way to do sensitivity analysis. This methodology clearly 
owns limitations i.e., results are strongly bounded to the 
values of fixed parameters, and different sets of values for 
the fixed parameters may entitle completely different 
results. 

Partial rank correlated coefficients (PRCC) [25] is a sta- 
tistical approach used to bypass the above mentioned 
limitations. It works by calculating the partial correlation 
on rank-transformed data between input (model para- 
meters) and output (entities behaviors). Such a technique 
does not depend on the values of fixed parameters and 
permits to vary all the parameters at the same time, 
allowing to study the influence of input parameters on 
the model outcomes. Nevertheless the methodology can 
be in principle easily applied and used with any kind of 
continuous or discrete model. 

The methodology we used to perform sensitivity analysis 
(LHS-PRCC) is briefly described as follows. The interested 
reader can found more information about the methodol- 
ogy in [26] . Parameters space is initially sampled using a 
Monte-Carlo technique. In this case we use a technique 
named Latin-Hypercube-Sampling (LHS) [27]. The techni- 
que divides the random parameter distributions into N 
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Table 2 Model parameters 


Param 


Description 


Value(estimate) 


Ref 




VC (VQ death rate 


ln(2) = 9 


In vivo 


a, 9 


VC killing rate by TC cells (TQ (killing) 


0.001 


Estimated 


a 17 


VC killing rate by AB (AB) molecules (killing) 


0.001 


Estimated 


q 


No. of cancer cells to inject at every vaccine administration 


50 


SimTriplex 


721 


released TAA (TAA) rate by killed VC 


3 


Estimated 


728 


released TAA rate by killed CC (CQ 


3 


Estimated 


to 


TAA natural degradation rate 


ln(2) = 9 


In vivo 


a2o 


Binding rate between TAA and APC cells 


0.0005 


Estimated 


a 27 


Binding rate between TAA and AB (IC formation) 


0.00001 


Estimated 


734 


plasma B cells (8) activation rate by TH cells (TH) 


0.05 


Estimated 


a 3 6 


B stimulation rate by IL2 (IL2) 


0.0035 


Estimated 




B duplication stimulation threshold due to IL2 


400 


Estimated 


to 


B cells natural death rate (half life) 


ln(2) = 15 


[32] 


740 


TH cells (TH) activation rate by APC (APQ cells 


0.15 


Estimated 


a 46 


TH cells stimulation rate by IL2 (IL2) (duplication) 


0.009 


Estimated 


Si 


duplication stimulation threshold due to IL2 


1000 


Estimated 


a 45 


TH cells cells stimulation rate by (//J 2) IL12 (duplication) 


0.009 


Estimated 




duplication stimulation threshold due to IL12 


1000 


Estimated 


to 


TH cells natural death rate (half life) 


ln(2) = 15 


Estimated 


751 


IL12 molecules release rate by VC 


10 


SimTriplex 


a 54 


absorbed IL12 rate by TH cells for mitotic signals 


0.00009 


Estimated 


a 59 


absorbed IL12 rate by TC cells for mitotic signals 


0.001 


Estimated 


to 


IL12 molecules natural degradation rate 


ln(2) = 9 


[33] 


764 


IL2 release rate by TH 


5 


Estimated 




absorbed IL2 rate by B cells for mitotic signals 


0.0001 


Estimated 


a 6g 


absorbed IL2 rate by TC cells for mitotic signals 


0.0001 


Estimated 


to 


IL2 molecules natural degradation rate 


ln(2) = 3 


Estimated 


773 


Released AB molecules rate by B cells 


3 


SimTriplex 


«78 


AB - CC binding rate 


0.0001 


Estimated 


a?i 


AB - VC binding rate 


0.001 


Estimated 


a?2 


AB - TAA binding rate (IC formation) 


= a?2 


_ 


to 


AB natural degradation rate 


ln(2) = 7 


[34] 




CC (CQ growth saturation threshold 


10 7 


Estimated 


/t 


CC duplication rate 


0.0226 


SimTriplex 


P 


No. of newborn CC due to transgenic nature of mice 


3 


SimTriplex 


ass 


CC death rate due to other IS entities 


0.0000001176 


Estimated 


a 89 


CC killing rate by TC cells 


0.00004 


Estimated 


as? 


CC killing rate by AB 


0.00004 


Estimated 


791 


TC cells activation rate by VC 


0.2 


Estimated 


a 96 


TC cells duplication rate due to IL2 


0.05 


Estimated 


596 


duplication stimulation threshold thanks to IL2 


400 


Estimated 


to 


TC cells natural death rate 


ln(2) = 21 


[35] 


702 


APC (APQ activation rate due to TAA (TAA) 


0.07 


Estimated 


to 


APC natural death rate 


ln(2)/15 


[36] 



(where N represents the chosen sample size) equal prob- 
ability intervals that are then sampled. The choice for N 
should be at least k + 1, where k is the number of para- 
meters varied, but usually much larger to ensure accuracy. 
In our trials we set N = 1000. 



After sampling an LHS matrix X of sampled parameters 
is built. Each row represents an unique set of variables 
for the model sampled without replacement. 

The model is then solved for each row of X, and the 
model output values are stored into an output matrix Y. 
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Each matrix is then rank-transformed (X R and Y R ). 
X and Y can be used to calculate the Pearson correlation 
coefficient. X R and Y R can be used to calculate the 
Spearman or rank correlation coefficient (RCC) and the 
partial rank correlation coefficient (PRCC). 

PRCC between an input parameter Xj e X R , j < k and 
output y e Y R is then computed by considering the resi- 
duals Xj — Xj and y — y where Xj and y are given by the 
following regression models: 



Xj = Co + 



k k 

Y] c p x p and y = bo + ^2 bpXp 
P=\,Pii p=i,prf 



Results and discussion 

The outcome of the in vivo experiment has been mainly 
represented by the mice tumor-free survivals, and the 
Kaplan Meier survival curves [4] for each mice group 
treated with a given vaccination protocol have been 
built accordingly (see Figure 3). 

One of the first problems in modeling the process was 
to determine how to translate the biological concept of 
death in mathematical/computational terms. When 
developing the SimTriplex model, it has been decided to 



stop the simulation and, therefore, to consider a mouse 
as dead if the total number of cancer cells reached 10 
cells. Over such a threshold the formation of carcinoma 
in situ can be considered an inevitable circumstance. 
Since in vivo experiments demonstrated that the vaccine 
progressively loses its efficacy when such an event 
occurs [21], this threshold represents a point of no 
return that halves between survival and death. 

Carcinoma in situ formation entitles a lot of different 
processes such as formation of physical barriers around 
the tumor mass and vascularization processes that are 
not described at this stage, even because this goes 
beyond the scope of the model. This means that both 
the ODE-based model and the SimTriplex models can- 
not be considered accurate in describing the in vivo 
experiment if the cancer cells threshold is overcome. 
Therefore the numerical simulations presented here 
refer to interactions where the number of cells do not 
go beyond this threshold. 

As previously stated, the success of failure of a treat- 
ment has been determined mainly by the survival rates 
of the mice involved in the experiment. Even if some 
measurements were made during the in vivo experi- 
ment, it was not possible to keep track of the time 
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Figure 3 Kaplan-Meier survival curves. Kaplan-Meier survival curves given by the in vivo experiment for the Untreated (red circles), Early 
(purple triangles) and Chronic (blue squares) vaccination protocols. 
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evolution of the involved entities. Such measurements 
are not possible in vivo experiments, or can be achieved 
just partially in vitro for multiple reasons, i.e. it is not 
possible to do the measure too frequently due to wet- 
lab requirements, it is not possible to take the measure 
at present time with current technology, or simply 
because the measure entitles the need to kill the host. 

One of SimTriplex main features is represented by the 
possibility of simulating different individuals. Tuning of 
free parameters has been executed in order to reproduce 
the same population survival curves for the vaccination 
protocols tested in vivo [5]. Moreover, during its tuning 
phase, SimTriplex entities behaviors have been accu- 
rately checked by biologists in order to verify that they 
were qualitatively in line with both biologists assump- 
tions and last immunological knowledge. The use of 
SimTriplex as a predictive tool, in conjunction with var- 
ious optimization techniques [28-30], to find better vac- 
cination protocols showed indeed that it represented a 
good approximation of the in vivo experiment [23], and 
therefore can be used to substitute missing in vivo data. 

Bearing all the above in mind, we initially checked 
that the mathematical model mice survivals for all the 
tested vaccine protocols were in tune with mean survi- 
vals showed in the in vivo experiment, obtaining a good 
agreement between the two experiments. For the miss- 
ing in vivo data, mainly represented by entities time- 
behaviors, we compared ODE behaviors obtained 
numerically with the ones obtained by SimTriplex, high- 
lighting similarities and differences. 

We would note here that, in order to compare the 
results, we looked in SimTriplex for "mean virtual 
mouse", i.e. a mouse whose death occurs near the middle 
of the Kaplan Meier curves for the tested protocols. 

The ODE model demonstrated able to reproduce the 
available in vivo experimental data, in particular the in 



silico mice survivals for all vaccine protocols tested were 
in good agreement with mean survivals showed in in 
vivo experiment. 

Since the biological behavior of the involved entities 
may change in a consistent manner even from mouse to 
mouse, we mainly focused in qualitatively analyzing can- 
cer cells behaviors and the response times of the princi- 
pal outcomes of immune response, i.e. antibodies and 
cytotoxic T cells behaviors for the Chronic, Early and 
Untreated protocols. 

In Figure 4 we compare the number of cancer cells 
(CC) behavior for the three protocols. As the Figure 4 
shows, there is a slight delay between SimTriplex and 
the ODE model curves for both Chronic and Early 
protocols, whereas the Untreated protocol exhibits 
negligible difference, since both models use the same 
parameters for the growth law. Such a delay remains 
in line with in vivo experiment expectations. Indeed 
the behaviors are qualitatively in agreement, suggesting 
that the cancer cells dynamics is well described by the 
ODE model. The Chronic protocol (see Figure 4, left 
panel) plot suggests that after an initial growth phase, 
cancer cells are kept under control from the immune 
system thanks to the repeated administration of the 
vaccine. 

The above behavior does not happen in the Early case 
(see Figure 4, center panel), where the vaccination pro- 
tocol is only able to delay the development of the can- 
cer, and the threshold on the number of cancer cells 
that entitles high risks of carcinoma in situ is reached at 
around at 44 weeks of age in SimTriplex, and at 47 
weeks in the ODE model. In in vivo experiments the 
middle of the Kaplan Meyer survival curve for the early 
protocol is reached approximately at 52 weeks of age 
[4], with carcinoma in situ formation between 5 to 9 
weeks earlier. 
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The number of cytotoxic T cells (TC) behavior is 
shown in Figure 5. The untreated plot (see Figure 5, 
right panel) is flat for both the models since in absence 
of vaccination there is no cytotoxic T activation. Even in 
this case we observe that the ODE model plots for the 
Chronic and Early vaccine protocols are a little bit 
delayed with respect to the SimTriplex plots, even if 
they remain in the expected range of the in vivo experi- 
ment. This could partially justify the delays observed in 
the cancer cells plots for the ODE plots (see Figure 4, 
left and center panels). 

Moreover the cytotoxic T cells peaks observed in Sim- 
Triplex for both the Chronic and Early protocols (see 
Figure 5, left and center panels) are a lot higher than 
those showed by the ODE model. In in vivo experiments 
it was observed that antibodies covered a major role in 
eradicating the tumor, whereas cytotoxic activity was 
estimated to be of secondary importance [23]. So from 
this point of view the ODE model may indeed be more 
accurate in describing this aspect of the immune 
response. 

Finally we analyze the number of antibodies {AB) 
behavior in Figure 6. The untreated plots (see Figure 6, 
right panel) are practically flat for both the simulations 
and show negligible difference. Only at the end of the 
experiment SimTriplex plots show the appearing of anti- 
bodies. This may be due to the presence in SimTriplex 
of Natural killer cells that are able to kill cancer cells 
which under-express the MHC molecules, giving rise to 
a week response at late stages. For the Chronic and 
Early protocols (Figure 6, left and center), the antibodies 
time-behavior is very similar, but with higher AB peaks 
in the ODE model than those observed in SimTriplex 
simulations. This can be seen as a consequence of the 
weaker cytotoxic immune response observed in the 
ODE model which requires, in accord with in vivo 



observations, a stronger humoral immune response in 
order to deplete the tumor. 

We used PRCC to analize the effects of the most 
important input parameters which influence more the 
behavior of Cancer Cells. We plotted for these entities 
the PRCCs over the entire time course of the experiment 
to how the parameters sensitivity varies as the process 
behavior advances. The analysis has been executed by 
supposing that the administration of the vaccine follows 
the Chronic protocol. In this way it is possible to study 
which mechanisms mainly drive the immune response 
against cancer cells and which parameters should be 
tweaked in vivo in order to obtain a strong immune 
response with the minimal effort. To this end we kept 
constant the parameter related to the quantity of injected 
vaccine cells (q) and the parameters related to the tumor 
growth (k, p, c max ). 

From the LHS-PRCC analysis we found that 15 para- 
meters that correlated significantly with the number of 
cancer cells. For some parameters a negative or positive 
correlation was somewhat expected, for example it is tri- 
vial to observe that the dead rate and the activation rate 
of APC (ft 0 and 702) positively and negatively correlate 
with cancer cells behavior, respectively (see Figure 7). 
The time correlation of some parameters indeed brings 
out some interesting findings we show as follows. 

The first interesting finding is related to the mechanisms 
driving the immune response against the tumor. The cyto- 
toxic immune response against the tumor is influenced by 
parameters such as 791 which represents the vaccine cells 
killing rate and a 89 which represents the cancer cells kill- 
ing rate by cytotoxic T cells. By taking a look at y 91 and 
a 89 (Figure 8) PRCC time plots it is possible to observe 
that the two parameters show a strong negative correlation 
just at the beginning of the experiment. The correlation 
becomes weaker and weaker as time goes on, becoming 
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Figure 6 Number of antibodies (AB) behaviors for the Chronic, Early and Untreated protocols Blue solid lines identify SimTriplex 
simulations, red dashed lines the ODE model numerical results. 



totally not significative starting from (around) day 160. 
The humoral immune response instead is driven by para- 
meters such as 7 73 which represents the antibodies release 
rate and a 8 7 which represents the cancer cells killing rate 
by antibodies. These parameters show a strong negative 
correlation which grows fast and lasts up to the end of the 
experiment (see Figure 9). This means that after an initial 
stage needed to start the immune response, variations of 
cytotoxic T cells related parameters do not influence can- 
cer cells behavior, thus suggesting how cytotoxic T cells 
are not fundamental for the complete eradication of the 
Tumor, which is instead strongly correlated with humoral 
immune response related parameters for all the time 
length of the experiment. This fact confirms the observa- 
tion made by Palladini et. Al. [23], where in vivo 



observations showed that the immune response against 
the tumor was mainly driven by antibodies. 

Another interesting finding can be derived by taking a 
look at 7 2 i and a 72 PRCC time plots (see Figure 10). The 
7 2 i parameter represents the antigens release rate by vac- 
cine cells whereas the a 72 parameter represents the rate 
of interaction between antibodies and antigens which 
brings to the formation of immune-complexes. The 721 
plot clearly shows a negative correlation with the number 
of cancer cells. However it is interesting to observe that 
this correlation becomes weaker at the end of the experi- 
ment when the vaccine protocol is used to maintain 
under control the number of cancer cells, thus suggesting 
that the role of antigens becomes less important at the 
end of the experiment. In addition the a 72 PRCC plot 
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Figure 7 u 0 and y 02 PRCC time plots. Partial Rank Correlated Coefficients are computed on the number of cancer cells (CQ, and are plotted 
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significant (p <0.01) are shown in gray. 



Bianca et al. BMC Bioinformatics 2012, 13(Suppl 17):S21 
http://www.biomedcentral.com/1471-2105/13/S17/S21 



Page 1 3 of 1 5 




at 

nt 

M 
o» 

! o 

-0» 

■«• 



<n) I'RCX" timt-plul fur ->,„ (b) l»RCC linn-plot (<w 
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shows a positive correlation with the number of cancer 
cells, suggesting that a higher interaction rate between 
antibodies and antigens negatively influences the effects 
of the treatment. This becomes more evident at the end 
of the experiment (just around day 300) when the last 4 
vaccination cycles are administered. Every vaccination 
cycle seems to reinforce the correlation between the 
input and the output parameters, affecting negatively the 
immune response. This can be explained by the fact an 
higher interaction rate between antibodies and antigens 
would entitle that more antibodies are recruited in bind- 
ing antigens, and then fewer antibodies (which as 



discussed earlier play a major role in the immune system 
response against the tumor) are employed in killing can- 
cer cells. When the number of cancer cells is kept under 
control, antigens may enter in competition with cancer 
cells and may negatively influencing the immune 
response. From the sensitivity analysis we can conclude 
that the antigenic administration should be then stronger 
during the first phases of the immune response, and then 
reduced once the humoral immune response is well 
established in order to reduce the risk that too many 
antibodies are involved in binding antigens instead can- 
cer cells. This model speculation is in line with the 
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Figure 9 y 73 and a 87 PRCC time plots. Partial Rank Correlated Coefficients are computed on the number of cancer cells (CQ, and are plotted 
over time (blue lines). PRCC plot of Dummy parameter (red lines) is presented for comparison. The plot portions where the correlation becomes 
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biological effect of high antigen stimulation that usually 
suppresses the immune response [31]. 

Conclusions 

The mathematical model proposed in this paper is based 
on nonlinear ordinary differential equations. The model 
simulates the competition between the immune system 
and the mammary carcinoma under the action of an 
external force field (the vaccine). Three different proto- 
cols of the vaccine have been taken into account: 
Untreated, Early, and Chronic. The biological role of vac- 
cine cells, cancer cells, tumor associated antigens, plasma 
B cells, thymus cytotoxic lymphocytes, thymus helper 
lymphocytes, antibodies, interleukins 2 and 12, and anti- 
gen presenting cells has been taken into account. 

Numerical simulations of the model have been per- 
formed for different vaccination protocols and results 
were compared with a previously developed multi-agent 
model, called SimTriplex. For the tested vaccination pro- 
tocols, the ODE-based model is able to qualitatively repro- 
duce the time evolution not only for the number of cancer 
cells, but also for antibodies and cytotoxic T cells, main 
outcomes of humoral and cell mediated immune 
responses. From a quantitative point of view the mathema- 
tical model showed, respectively, a weaker and a stronger 
immune response of cytotoxic T cells and antibodies with 
respect to the SimTripex model, showing indeed better 
agreement with the in vivo observations and speculations. 

The sensitivity analysis gave two major results. First it 
confirmed the major role of humoral immune response 
also observed in in vivo experiments [23], then showed 
that during later stages of the experiment antigens loose 
their role of activating the immune response and in 



some cases may negatively influence the immune 
response. It is then possible to conclude that a reduction 
of the intensity of vaccine administrations in later stages, 
when the immune response is already set, is advisable. 
This has been also highlighted in [8], where an ABM 
model developed to illustrate the effects of the same 
vaccine in cancer immunotherapy, suggested to apply 
the golden standard vaccination procedure (initial boost 
followed by sparse recalls) also to cancer vaccines. 

These results are certainly useful to research activity in 
immunology addressed to improve the efficacy of the treat- 
ment and to modulate the activation of the immune system 
in order to prevent side effects such as autoimmune dis- 
eases. Of course, different choices of initial conditions and 
of the parameters may modify the competition dynamics. 

We plan to investigate the optimal protocol using 
mathematical tecniques which are currently under inves- 
tigation. Results will be pubblished in due course. 
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